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I . Purpose 


The purpose of this progress report is to summarize the 
research progress for grant NAG-1-1350 entitled "Finite 
Difference Time Domain Modeling of Steady State Scattering from 
Jet Engines with Moving Turbine Blades" for the time period 
9/25/91-3/25/92. 


II. Introduction 


Steady state scattering, especially determination of Radar 
Cross Section, from jet engines with moving turbine blades is a 
difficult problem. Doppler frequency components are to be 
expected in the frequency domain scattering response and could be 
used to detect the rotational speed of the engine and hence 
possibly the type of engine. This is of great significance to 
detection and target identification for current and future 
defense/navigational radar systems. 

The approach chosen to model steady state scattering from 
jet engines with moving turbine blades is based upon the Finite 
Difference Time Domain (FDTD) method. The FDTD method is a 
numerical electromagnetic program based upon the direct solution 
in the time domain of Maxwell's time dependent curl equations 
throughout a volume. One of the strengths of this method is the 
ability to model objects with complicated shape and/or material 
composition. General time domain functions may be used as source 
excitations. For example: a plane wave excitation may be 
specified as a pulse containing many frequencies and at any 
incidence angle to the scatterer. A best fit to the scatterer is 
accomplished using cubical cells in the standard cartesian 
implementation of the FDTD method. The material composition of 
the scatterer is determined by specifying its electrical 
properties at each cell on the scatterer. Thus, the FDTD method 
is a suitable choice for problems with complex geometries 
evaluated at multiple frequencies. It is assumed that the reader 
is familiar with the FDTD method. References [l]-[3] provide 
general background information on the method. 

The challenges are many in this analysis. Upon first glance 
the most difficult task appears to be the integration of moving 
objects into a stationary analysis. However, other factors must 
also be considered. Namely, the jet engine appears to be a 
resonant cavity with a high Q. This means that the 
electromagnetic energy that penetrates the inlet stays inside for 
a very long time. This in turn requires a much longer 
computation time in order to reach the steady state field value. 
Another obstacle must be overcome in addition to the two 
impediments mentioned already: the size of the scatterer. The 
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issue here is not with the volume the scatterer occupies within 
the problem space but with the shape and the dissimilar 
proportions of the components of the scatterer. These three 
issues are discussed in the following sections. 


III. Engine Geometry 


A brief review of the FDTD method has already been 
conducted. This section is concerned with the scattering object 
that is constructed in the problem space. The scatterer is 
located in the volume by discretizing the object into cubical 
cells and assigning electrical characteristics (permittivity, 
conductivity) to each individual cell in the x, y, and z 
directions. The scatterer is easily modified simply by changing 
the electrical characteristics at individual cell locations. The 
ability to change the scattering object (the volume) is one of 
the strengths of the FDTD method. 

One of the difficulties inherent to FDTD is the problem 
space size. The issue here is twofold. The scatterer is large 
and the actual volume the scatterer occupies within the problem 
space is appreciable. Additionally, the shape and dissimilar 
proportions of the components of the scatterer are of concern. 

And still other issues must be considered based on the 
characteristics of scatterer, both electrical and geometrical. 

First let us examine the size of the problem space. In 
order to size the problem one must consider the machine on which 
the problem may be evaluated. The computer on which this 
calculation is based is a Silicon Graphics workstation with 64 MB 
of RAM. The standard Yee unit cell (Figure 1) is used and 
requires six real field component variables as well as three 
integer* 1 material property variables for characterization of 
perfect electrical conductor, dielectric and lossy dielectric 
structures. This yields a total of nine variables per cell. 

Each field variable requires 4 bytes of memory and each material 
property variable requires 1 byte of memory so 27 bytes are 
required to store variables for each cell. The limit of RAM 
accessible for a computation on this workstation is approximately 
60 MB since 4 MB must be reserved for system utilization. 

In order to maximize utilization of the problem space volume 
the geometry of the scatterer must be considered. Since aircraft 
engines are longer than they are wide, the engine inlet is chosen 
to be four times as long as it is wide, Figure 2. The cell 
volume of the scatterer may be expressed as x by x by 4x whxch 
equals 4X 3 . The variable x represents the total number of cells 
across the short dimension. Now let us determine the number of 
cells along the short dimension of the problem space x: 
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Volume - 


Available RAM 
storage per ceTl 


A 3 60 * 2 20 bytes 

qx ' 27 bytes 


4x 3 = 2.33 * 10 6 


=> x = 83.5 cells . 


Assuming that a 10 cell border must surround the scattering 
object, the available cells for modeling are 73 X 73 X 292 - 
1,556,068. For a small engine with a diameter of 1.0 m the cell 

size Ax is 

Av Engine diameter 
~lno3eTTng~ceTTs' 

. 1.0m 

Ax ' 73 cells 

=> Ax = 1.37 cm. 


With Ax = 1.37 cm and a resolution of 10 cells per wavelength the 
maximum frequency allowable is 2.19 GHz. 

Not only is the scatterer very large, but it also contains 
fine components , see Figure 3 . The moving parts of the engine 
are of interest in this analysis and the smallest components, 
such as the hoses and wires, will not be considered at this time. 
The dilemma is to accurately render the turbine blades when 
restricted to a minimum cell size by the computational machine 
limits . 

Other difficulties occur when the geometrical configuration 
of the scatterer is a semi-closed shape. Contrary to intuition, 
electromagnetic energy penetrates a semi-closed structure and, 
depending on the loss characteristics of the material, remains 
inside for long periods of time. Once the energy is inside of 
the semi-closed structure the frequencies of the penetrating 
signal are re-tuned to the resonant frequencies of the structure 
[4], Thus the semi-closed structure is resonant at certain 
frequencies dictated by the geometrical configuration of the 
structure. The quality factor, Q, relays the resonant 
information of a structure. High Q values indicate extremely 
resonant cavities. The high Q of the engine means that the 
computation must be run for a greater number of time steps if the 
conventional post-processing approach is used on the time domain 
results . 
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The FDTD method has been applied to resonant cavities [5]. 
Based on this analysis and the geometrical configuration and 
electrical properties of the engine, it is likely that the engine 
is a highly resonant structure. The engine appears to be a semr- 
closed structure due to the overlapping blades of the compressor 
at the front and the turbine at the back of the engine. This 
combined with the low loss material (metal) of the engine 
components cause the scatterer to be a resonant structure . 
Therefore electromagnetic energy within the scatterer lingers tor 
a long time. Hence, a significant number of time steps will be 
required. Using the Ax of 1.37 cm calculated previously the time 
step At is 


At = 


Ax 




At = 


1.37 cm 


y/J (3.0 * 10 8 m/s) 


=> At = 26.4 ps. 

Estimating that most of the energy has dissipated by 1 ns, the 
number of time steps required is 

Cavity decay time 
TIME STEPS = — l — 


TIME STEPS = 


1.0 * 10- 06 


26.4 * 10- 12 
TIME STEPS = 37,928. 


The cavity decay time may be smaller than this when the effects 
of loss in the engine, e.g. carbon from burning fuel, are 
included. 

Having discussed the size of the scatterer and its 
predisposition to retain energy, a first attempt at depicting the 
engine geometry is shown in Figure 4. Figures 5, 6 and 7 show 
the front view, side view and the centerline slice of the engine. 
This is a zeroth order model and serves as a good starting point 
for determining the resonance of the engine and the effectiveness 
of the Motion FDTD computation in evaluating the rotating parts 
of the engine. The zeroth order model has been verified using a 
three dimensional FORTRAN graphics program on the Silicon 
Graphics workstation. 

The zeroth order model is constructed in three main parts . 
The outside of the engine is called the inlet and is circular 
cylindrical in shape and hollow inside. The inlet is the 
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stationary part of the engine. The moving components are the 
rotor shaft and the blades. The shaft is modeled as a solid 
circular cylinder which extends from the front to the rear of the 
inlet. The blades form a planar fan in the zeroth order model. 
The fans are identical and one represents the compressor at the 
front of the engine and the other represents the turbine at the 
rear. The two-dimensional fans lie in the plane normal to the 
engine axis and intersect the axis at the front and rear of the 

inlet. The zeroth order model is automated and the user 

specifies the number of blades in the fan unit, the thickness of 
the blades, the radius of the rotor shaft, the radius of the 

blades, and the radius of the engine inlet. 

The number of blades specified for this model must be even. 
The blades are distributed symmetrically about the engine axis. 
The base of the blade starts on the shaft at r SHAFT and the tip of 
the blade extends to r BLADE . The width of the blade refers to the 
area between the two radii that may be metal. The width or tne 
blade controls the openness of the engine. By specifying the 
number of blades and the width of metal to air of the blade arc, 
the model serves to approximate the overlapping blades of a 
turbine/compressor assembly. In this way the ratio of air to 
metal in the planar fan remains constant but the occurrence of 
the air to metal pattern increases and the arc length subtended 
by each decreases as the number of blades increases . 

The number of blades, NB, is specified and a new blade 
starts every 0 incr where incr stands for increments and 

0 - 2jt 
incr- NB * 

The ratio of metal to air, BLDWD, is specified as a real number 
between 0.0 and 1.0. If a zero is specified the fan structure is 
not constructed, i.e. the blades have zero thickness so the fan 
is air. When a one is specified the fan structure becomes an 
solid metal disk (like a coin), i.e. the blades completely fill 
the planar fan with metal. If 0.5 is specified this means that 
half of the blade is metal and half is air. Thus for the first 
blade the metal starts at 0.0 degrees and continues until A0 
where 

A0 = BLDWD * & incr • 


The remainder of the angle from A0 to 0 incr is air. Figure 8 
illustrates these variables for one blade of the planar fan. 
Figure 5 depicts the blade structure when eight blades are 
specified and BLDWD is 2/3. 

Using the zeroth order model of the engine, we plan to 
investigate the Q of the cavity. The next model of the engine 
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will be based on the zeroth order model. The rotor and inlet 
will be the same and the fan will be extended to three 
dimensions. The decay rates of the fields in the zeroth order 
and first order model will be compared. This will be useful in 
determining whether a three dimensional Motion FDTD algorithm is 
in order or if a two dimensional Motion FDTD algorithm may be 
used. Note that if the fields in both models behave similarly 
then the zeroth order model of the turbine/compressor assembly 
might suffice. However if the scatter from the first order model 
using a three dimensional fan is markedly different from t e 
scatter from the zeroth order model, the next order model might 
have to be used. The effect of the three dimensional fan on the 
cavity fields of the engine will be investigated in the next time 

period. 


IV. Problem Space Termination 


It is evident from the discussion above that any 
computational savings is welcome. With this in mind, consider 
that additional cells must be added to the exterior of any 
scattering object evaluated using FDTD in order to lower the 
erroneous reflection from the problem space termination. A e 
very least a ten cell free space border must be attached to every 
dimension of the scatter when applying the second order absorbing 
boundary conditions developed by Mur [6]. However, when the 
scattered fields are highly reactive or when the scattered fields 
are not normally incident to the absorbing cells, even a ten cell 
border does not prevent large reflections from the problem space 
termination. 

An improvement to the absorbing boundary conditions 
developed by Mur has been tested. The improved boundary 
conditions are called Super Outer Radiating Boundary Conditions, 
Super ORBC , or SORBC [7]. The SORBC minimize the non-normal 
incidence and reactive errors mentioned above. 

Determination of the electric and magnetic field at the 
termination of the computational grid using the SORBC involves 
several steps. First, the incidence angle of the scattered field 
at each border cell is determined from the Poynting vector. 

Next, electric and magnetic fields are estimated as follows. The 
first order Mur equations are used to predict edge electric field 
components from the inner electric field components. The same 
equations are then used to predict an edge magnetic fxeld 
component. The estimated magnetic field component is used by the 
finite difference equations to estimate yet another electric 
field value. This results in two different estimates of the 
electric field for each border cell. Finally, both estimated 
electric field values are utilized using a weighted average. The 
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two estimated electric fields are added using a weighted average 
which depends on the arrival angle of the Poynting vector. 

The SORBC provides more accurate results than the second 
order Mur equations because the incidence angle is included in 
the calculation. The SORBC have been verified and will be 
described in more detail in the next progress report. 


V. Moving Objects 


One of the difficulties in this problem is because of the 
movinq objects that are present in the analysis. An approach 
needs to be developed that properly treats the physical phenomena 
involved with scattering from moving objects and this approach 
needs to be adaptable to the FDTD method. Another difficulty 
arises in incorporating any approach into the FDTD method because 
of the moving objects. In the FDTD method, the object and 
surrounding computational space are sampled at discrete points in 
space and electric and magnetic fields are interleaved on this 
spatial grid. Once an object begins to move, it is no longer 
"attached" to the grid and becomes more difficult to analyze 
using FDTD. A discussion of the two main approaches investigated 
follows . 


A. Single Grid Approaches 


The first approach tried for simulating moving objects in 
FDTD was to allow the object to move one cell every N time steps. 

N is calculated from the duration of one time step to give the 

desired speed: eg. N = (Ax/v x )/At, where v x is the x-component 

of velocity. This is referred to as the " jumping— cell ^algorithm 
because the moving object appears to suddenly pass or jump from 

one cell to the next every N time steps. The fields and 

constitutive parameters associated with a moving cell are copied 
to the new location so there is a crude approximation of charge 
flow along with motion of the object boundary. In spite of 
approximating the charge flow at the object boundary, this 
technique is still very non-physical. One reason is the 
relatively large discontinuity that occurs between the two 
adjacent time steps every time the object moves. Another reason 
is the failure to model the instantaneous field changes in front 
of the object. Figure 9 shows the results of this simulation on 
the E-field backscattered from a perfectly conducting slab. In 
this figure, the scattered field from a stationary slab is 
compared with that from the same object in "motion." The slab 
velocity in the moving case is c/100 m/s, which should give a 
frequency shift of about 10 MHz. No shift is noticed, though 
sidebands are produced at harmonics of the cell moving rate. 
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The second approach, referred to here as the "growing-cell 
algorithm, is a significant refinement of the above procedure. 
The "qrowing-cell" approach allows the object to move a fraction 
of a cell length every time step. Consider the free space cell 
directly in front of a moving conductor, where the conductor 
velocity is normal to its surface. This is shown in Figure 10. 
The integral form of Maxwell's electric curl equation expressed 
around the perimeter of this cell is 

f d = - -(£)// »• d_s = + J > 1 dl > • 

In finite difference form this gives (for velocity in the x- 
direction) 


Tjfi + X/ 2 1 / 2 


At (A E y _ A E z \ 

p. {~£z ~Sy ) 


^3" < 


H 


n- 1/2 


»»Jl + l/2 rjM— 1/2 

ti y — tly 


At [AF z 
| x Ax 



+ -^H n y - 1/2 • 


Results for this approach are given in Figure 11 for the same 
problem geometry and slab velocity as in Figure 9. The moving 
slab case shows the correct Doppler frequency shift, although 
sidebands are still produced and the Doppler amplitude, which 
should have increased less than a dB, is incorrect. 


B. Dual Grid Approach 


The Dual Grid approach involves a Lorentz type of 
transformation of the field components and interpolation into the 

FDTD grid. 

A recent paper [8] describes electromagnetic scattering ftom 
moving surfaces in one and two dimensions using the FDTD method. 
In that paper, the motion of the surfaces were tracked in the 
FDTD code and the relativistic boundary conditions were 
implemented at each surface and the total fields were set to 0 
inside each conductor. The results in the paper agreed 
reasonably well with published results, but the limitations are 
that only perfect conductors were considered and the analysis was 
not generalized. The paper also mentions the approach of using a 
Lorentz transformation of the incident field to the moving system 
and then solving the scattering problem in the moving system 
where the geometry is stationary. The final step involves a 
Lorentz transformation back to the inertial frame of reference. 
This is the basic approach that is developed for the present 
work. 
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Since jet engine turbine blades exhibit strictly rotational 
motion, our approach is based upon a transformation to and from a 
rotating system of coordinates. Many articles have been 
published regarding Maxwell's equations in rotating systems and 
the solution of scattering problems with rotating objects [9]- 
[191. Maxwell's equations are invariant under any coordinate 
transformation, but the constitutive equations are modified. Let 
us denote the laboratory frame of reference as F, the rotating 
frame of reference R and quantities in frame F will be unprimed 
and those in frame R will be primed. The rotating frame R has a 
constant angular velocity Q (see Figure 12). In the laboratory 
frame, Maxwell's curl equations are [16] 


V* E = 
Vx H = 


dB 

TE 


dD 

Tt 


+ j 


( 1 ) 


and the constitutive relations become modified as 

B = p [h- vx d] + _L V* E 


1 vx H 


D = E [fi + vx fl] - — _ 


J = ya |E + v x B - * -E 


( 2 ) 


+ pv 


where 

v = Q x r , y = = — i P = — ( 3 ) 


c is the speed of light in free space and p is charge density. 

The transformation to rotating coordinates is given by 

x' = xcos(Qt ) + ysin(Qt ) 
y' = -xsin(Qt ) + ycos(Qt ) 

z' = z ' 

t> = t 

(J)' = <j) - Qt 


and the inverse transformation is obtained from (4) by 
interchanging the primed and unprimed variables and replacing Q 
by -Q. Maxwell's curl equations remain invariant in rotating 
coordinates and are given by 
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v 7 x F = - — 

at' 

V'x?=^J 


(5) 


with the constitutive relations defined by 

B 7 = utf 7 + J_ vx E 7 
c 2 

D 7 = eE 7 - J_v* H 7 
c 2 

j 7 = oE 7 . 

The transformation of the field components is given by 

E* = E + vx B_ 

H 7 = H - vx D 
B 7 = B 
D 7 = D . 


( 6 ) 


(7) 


In this approach, two FDTD grids are used in the analysis. 
One grid is a stationary grid and is located in the laboratory 
(or inertial reference) frame and the other is a grid that is 
"attached" to the rotating object (see Figure 13). The 

stationary grid will be a minimum of cells larger on a side 
than the rotating grid so that all field components contained in 
the rotating grid will remain inside the stationary grid as the 
blades rotate. Thus, the object remains stationary with respect 
to R, and the FDTD grid points are well defined within the 
rotating frame. However, as the object rotates, the grid points 
in R will no longer coincide with the grid points in the 
laboratory frame F. The procedure to obtain the FDTD field 
components in frame F will involve interpolation of the field 
components in frame R and then a transformation back to frame F. 
The full procedure is outlined as follows: 

1. ) Transform the incident field to rotating coordinates using 

(4). 

2. ) Compute field values in the rotating frame using (5) and 

( 6 ). 

3. ) Apply boundary conditions in the rotating space R. 

4. ) Interpolate fields from the rotated space R spatial 

locations in R corresponding to FDTD field locations in the 
laboratory frame F. 

5. ) Transform interpolated fields from step 4 to the laboratory 

frame I using ( 7 ) . 

6. ) Apply stationary FDTD equations to field components in the 
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laboratory frame F that were not interpolated from the 

rotating space R. 

7.) Apply boundary conditions to the stationary grid. 

Step 1 of this procedure has been successfully demonstrated 
using a two dimensional TM code containing only the FDTD incident 
field. Two grids were defined in the normal FDTD manner and the 
incident field was analytically specified as 

E*(t) = r ( t) sin(w 0 t ) 

where r(t) is a ramping function that transitions from 0 to 1 
over several cycles of the sinusoid. For a test problem a 
sinusoidal frequency of 20 Hz was chosen so that Doppler 
frequency components would be very apparent. Of course, higher 
frequency incident fields and wideband pulses will be used in the 
full analysis. A point within the FDTD grid was chosen and the 
incident field was computed both in the stationary grid and the 
rotating grid. For 0=0, both incident fields were identical as 
expected. According to [17], the incident field in rotating 
coordinates should exhibit an infinite number of sidebands 
corresponding to 

(j) w = o> 0 + mCi ( ® ) 

for nonvanishing Q. This behavior does not become readily 
apparent in the frequency response until Q becomes on the order 
of to 0 . Figures 14 — 15 show the time and frequency domain response 
of the incident field in stationary and rotating coordinates for 
Q=2 ji rad/sec corresponding to Q/co 0 =0.05. Note that the time 
domain rotating incident field is almost identical to the 
stationary incident field except for a slight time shift 
(corresponding to a phase shift). The frequency response of the 
rotating incident field is slightly shifted from the stationary 
response due to the smearing of the sideband components that are 
spaced only 1 Hz apart. Figures 16-17 show the time and 
frequency domain incident field for Q=20 ji rad/sec corresponding 
to Q/to 0 =0 . 5 . Note that the rotating time domain field is now 
markedly different than the stationary field and Figure 17 shows 
the frequency response indicating the sidebands spaced 10 Hz 
apart according to (8). These effects were also seen for larger 
ratios of Q/to 0 . 

The dual grid approach shows the most promise at this point 
for a computational approach to modelling rotating objects. The 
incident field in the rotating system exhibits sidebands at 
multiples of the rotational speed according to theory. The next 
phase of development is to test this approach on a two- 
dimensional problem similar to one analyzed in [19] and compare 
with another computational or analytic solution. If the approach 
works on a two-dimensional test problem, then it can be 
generalized to three dimensions to analyze an engine geometry 
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complete with canted rotor and turbine blades . An improvement to 
the brute-force dual grid approach which will also be 
investigated is to apply the dual grid technique on a local basis 
for only those portions of the engine that are rotating. This 
improvement would make the dual grid approach more efficient 
since two three-dimensional grids would not be used thereby 
decreasing memory and CPU requirements. 

VI . Summary 


Referring to the Statement of Work for "Finite Difference 
Time Domain Modeling of Steady State Scattering from Jet Engines 
with Moving Turbine Blades," many of the preliminary tasks have 
been completed or are near completion. The status of and a brief 
summary of the tasks attempted during this time period follow. 

Task 1: Automate Target Geometry Input 

An automated engine geometry computer code was developed. 
This model is called the zeroth order model because the 
turbine/compressor fan is confined to a plane. The user 
initializes several engine dimensions and the material 
properties at each FDTD cell are automatically assigned. 

Task 2: Graphical Verification of Target Geometry 

The graphical verification of the engine geometry was 
developed in conjunction with the zeroth order engine model. 
A three dimensional graphics program for visualization of 
the three dimensional engine geometry was completed. This 
program is written in FORTRAN and utilizes the graphics 
library on a Silicon Graphics workstation. 

Task 3: Graphical Display of Electromagnetic Fields 

A program to show the electromagnetic fields in the problem 
space is near completion. 

The absorbing boundary conditions at the edge of the 
computational grid have been improved. The modified grid 
termination equations, the SORBC, reduce reflection errors 
attributed to non-normal incidence angles and highly reactive 
fields . 

Several approaches to modeling the rotating objects in the 
jet engine using FDTD were considered. Further investigation of 
both the Single Grid Approach and the Dual Grid Approach is 
needed . 
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VIII. Figure Titles 

1. ) Standard three-dimensional Yee cell showing placement of 

electric and magnetic fields. 

2. ) View of engine inlet showing diameter d, length 4d and 

problem space dimensions x and 4x. 

3. ) View of jet engine showing various components. 

4. ) View of jet engine inlet showing zeroth order geometry 

model . 

5. ) Front view of engine geometry for an eight blade example 

with variable BLDWD=2/3 and showing various radii. 

6. ) Side view of engine geometry for the eight blade example of 

Figure 5 . 

7. ) Centerline slice of engine geometry for the eight blade 

example of Figure 5 . 

8. ) Zeroth order blade model showing variables 9i nC rf ^0, r shatt 

and r blade . 

9. ) Scattered electric field versus frequency for moving slab 

using "jumping-cell" algorithm with v=c/100 m/s. 

10. ) Cell in free space immediately in front of moving conductor 

showing integration loop used to derive FDTD H-field update 
equation for "growing-cell" algorithm. 

11. ) Scattered electric field versus frequency for moving slab 

using "growing— cell" algorithm with v=c/100 m/s. 

12. ) Inertial reference frame F and rotating frame R. R rotates 

with constant angular speed Q with respect to F. 

13. ) Two-dimensional stationary grid (unprimed quantities) and 

two-dimensional rotating grid (primed quantities) used for 
dual grid approach. 

14. ) Time domain incident field in stationary and rotating 

coordinates with co 0 =40ji rad/sec, Q=2it rad/sec and Q/to 0 =0.05. 

15. ) Frequency response of incident field in stationary and 

rotating coordinates with to 0 =40ji rad/sec, C2=2 ji rad/sec and 
Q/oj 0 =0 . 05 . 

16. ) Time domain incident field in stationary and rotating 

coordinates with (o 0 =40ji rad/sec, Q=20 ji rad/sec and Q/u) 0 — 0.5. 
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17*) Frequency response of incident field in stationary and 

rotating coordinates with a) 0 =40ji rad/sec, Q=20 jc rad/sec and 
Q/oj 0 =0 . 05 . 
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